%%%Main File for BCDR Simulations - gains from industrial policy and trade
%%%policy

%This file computes the gains from a given trade and industrial policy
%using the SOE assumption, with no intermediates

%CES Version

clear all;

%digits(64);

K=32; %number of sectors
N=62; %includes ROW
K_man=15;    %number of manufacturing sectors

%%%Load data
%%%Data is organized as follows: first column origin country codes, second column
%%%sector codes, third column sectoral value added, fourth column aggregate origin trade balance,
%%%fifth column origin-sector expenditurs, 6th column starts the bilateral
%%%trade flows

%Note: data is sorted by sector, then origin.

A=importdata('simulation_ge_data.csv');
B=importdata('sim_params_no_int.csv'); %these are the estimated elasticities

%Data vectors
X_ink = A(:,6:5+N); %Matrix of trade flows, (NxK)xN
D_i=A(1:N,4); %vector of trade balances, Nx1
V_ik=A(:,3); %vector of sectoral value added, (NxK)x1
expenditures_ik=A(:,5); %vector of expenditures, (NxK)x1
countryid=A(1:N,1); %vector of country ids, Nx1

%Parameter vectors
theta_nm=4.32; %non-manufacturing sector TE set to median of manufacturing
h_theta = theta_nm*[ones(2,1);zeros(K_man,1);ones(K-2-K_man,1)]+[zeros(2,1);B(:,2);zeros(K-2-K_man,1)]; %trade elasticities
h_gamma = [zeros(2,1);B(:,1);zeros(K-2-K_man,1)]; %Heterogeneous scale elasticities, non-manufacturing set to zero
gammaXtheta=h_theta.*h_gamma;
rho=.87; %Upper tier elasticity of substitution

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reshaping data vectors

%Reshape trade matrix to be NxNxK (3 dimensional array)
%rows are exporters, columns are importers, 3rd dimension is sectors
X=zeros(N,N,K);
for k=1:K
    X(:,:,k) = X_ink(1+(k-1)*N:k*N,:);
end
X_ink=X;

%Other data matrices
V_ik = reshape(V_ik, [N 1 K]);
V_i=sum(V_ik,3);
D_i =reshape(D_i, [N 1 1]);
expenditures_ik = reshape(expenditures_ik,[N,1,K]);
expenditures_nk=permute(expenditures_ik,[2 1 3]);

%Create trade shares, expenditure shares
expenditures_nnk=repmat(expenditures_nk,[N 1 1]);
lambda_ink = X_ink./expenditures_nnk; %trade shares

expenditures_n=sum(expenditures_nk,3);
beta_nk=expenditures_nk./repmat(expenditures_n,[1 1 K]);

%Create sectoral domestic trade shares
dom_lambda_ik=zeros(N,K);
for i=1:N
    for k=1:k
        dom_lambda_ik(i,k) = lambda_ink(i,i,k);
    end
end
        
%Create total sectoral exports
exports_ik=zeros(N,K);
for i=1:N
    for k=1:k
        exports_ik(i,k) = V_ik(i,k)- X_ink(i,i,k);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from Industrial Policy

%Defining output vectors
gfop_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)
gftp_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)
gfip_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)

%Output vectors for Harberger analysis
VA_opt=zeros(K,N);  % value added share by sector at the optimum
lab_real_opt=zeros(K,N);  %percent change in labor allocated to sector k when going from optimum to original equilibrium
fs_opt=zeros(K,N);   %foreign expenditures on domestic goods, as a share of domestic value added, at the optimum
flab_real_opt=zeros(K,N);  % (negative) percent change in effective labor units supplied to foreigners when going from the optimum to original equilibrium

VA_initial=zeros(K,N);  %value added share by sector at the laisezz-faire equilibrium
fs_initial=zeros(K,N);  %foreign expenditures on domestic goods, as a share of domestic value added, at the initial equilibrium

lab_real_tp=zeros(K,N);  %percent change in labor allocated to sector k when going from trade policy equilibrium to original equilibrium
flab_real_tp=zeros(K,N);  % (negative) percent change in effective labor units supplied to foreigners when going from the trade policy equilibrium to original equilibrium

lab_real_ip=zeros(K,N);  %percent change in labor allocated to sector k when going from industrial policy equilibrium to original equilibrium
flab_real_ip=zeros(K,N);  % (negative) percent change in effective labor units supplied to foreigners when going from the industrial policy equilibrium to original equilibrium

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculating Gains from trade as well
gft_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculating some correlations
corr_se_va_avg=zeros(N,1); %correlation between SE and value added
corr_se_lambda_avg=zeros(N,1); %correlation between se and import shares of expenditure
corr_se_exports_avg=zeros(N,1); %correlation between se and export shares of value added

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from trade and industrial policy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Het TE, Het SE (service 0 SE)
theta=h_theta; 
gamma = h_gamma;
counter=0;

for i=1:N
    %Defining data vectors and parameters
    V = V_i(i); %total value added
    X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
    D= D_i(i)/V; %Aggregate trade balance (divided by value added)
    beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
    lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
    V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
    
    VA_initial(:,i)=V_k; 
    fs_initial(:,i)=X;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Optimal trade and industrial policy
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gfop_hte_0(i,1)=model_output.welfare;
    
    %Saving reallocation terms for Harberger analysis
    
    VA_opt(:,i)=model_output.VA_share;  % value added share by sector at the optimum
    lab_real_opt(:,i)= model_output.lab_real;  %percent change in labor allocated to sector k when going from optimum to original equilibrium
    fs_opt(:,i)= model_output.fs_share;   %foreign expenditures on domestic goods, as a share of domestic value added
    flab_real_opt(:,i)=model_output.fl_reallocation;% (negative) percent change in effective labor units supplied to foreigners when going from the optimum to original equilibrium
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Trade policy only
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=zeros(K,1);
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    VA_TP=model_output.VA_share;
    gftp_hte_0(i,1)=model_output.welfare;
    lab_real_tp(:,i)=model_output.lab_real;  %percent change in labor allocated to sector k when going from trade policy equilibrium to original equilibrium
    flab_real_tp(:,i)=model_output.fl_reallocation;  % (negative) percent change in effective labor units supplied to foreigners when going from the trade policy equilibrium to original equilibrium
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Industrial policy only
    %Specifying taxes and subsidies
    t_k=zeros(K,1);
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gfip_hte_0(i,1)=real(model_output.welfare);
    lab_real_ip(:,i)=model_output.lab_real;  %percent change in labor allocated to sector k when going from industrial policy equilibrium to original equilibrium
    flab_real_ip(:,i)=model_output.fl_reallocation;  % (negative) percent change in effective labor units supplied to foreigners when going from the industrial policy equilibrium to original equilibrium
    
    %Total reallocation for gains from industrial policy
    tot_reallocation (i,1) =sum(abs(VA_opt(:,i)-VA_TP)); 
    
    %Gains from trade
    model_output=gft_ces(V_k,lambda,beta,gamma,theta,rho,K);
    gft_hte_0(i,1)=model_output.welfare;
    
    %Some correlations
    z= corrcoef(V_k,gamma);
    corr_se_va_0(i,1)=z(1,2);
    z= corrcoef((beta.*(V+D*V)).*(1-lambda),gamma);
    corr_se_lambda_0(i,1)=z(1,2);
    z= corrcoef(X,gamma);
    corr_se_exports_0(i,1)=z(1,2);
    %Now just the product of SE and value added
    prod_se_va(i,1)=V_k'*gamma;
    
    counter=counter+1
end

%%%%%%%%%%%
%Computing gains from industrial policy
gfoip_hte_0=gfop_hte_0-gftp_hte_0;

%Saving the gains
gains_hte_0=[countryid gfop_hte_0 gftp_hte_0 gfip_hte_0 gfoip_hte_0 gft_hte_0  corr_se_va_0 prod_se_va tot_reallocation  corr_se_lambda_0 corr_se_exports_0];
csvwrite('gains_soe_hte_ces_0.csv',gains_hte_0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Harberger formula with open economy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gamma=h_gamma;

%Open Economy - NEW FORMULA
% Create change in sector scale and change in effective labor devoted
%to exports when going from equilibrium with optimal policy to
%laissez-faire equilibruim.  Use these to compute the formula.

%Note: we have a few zero export sectors, for which the formula delivers an NAN
%value for foreign labor reallocation.  Need to replace with zero (which is
%appropriate)
flab_real_opt(isnan(flab_real_opt))=0;
flab_real_tp(isnan(flab_real_tp))=0;
flab_real_ip(isnan(flab_real_ip))=0;

lab_real_opt_tp=zeros(K,N);
for i=1:N
    reallocation_opt =1./(1-lab_real_opt(:,i))-1;   %l_hat -1 when going from initial to optimum
    fl_reallocation_opt =1./(1-flab_real_opt(:,i))-1;   % percent change in effective labor units supplied to foreigners when going from the original equilibrium to optimum
    dw_harberger_opt(i,1) =(1/2)*sum(VA_initial(:,i).*reallocation_opt.*gamma+fs_initial(:,i).*fl_reallocation_opt.*-1./(theta+1));% Harberger gains from optimal policy
    dw_harberger_opt_trade(i,1)=(1/2)*sum(fs_initial(:,i).*fl_reallocation_opt.*-1./(theta+1)); %Optimal gains due to trade policy
    dw_harberger_opt_industrial(i,1)=dw_harberger_opt(i,1)-dw_harberger_opt_trade(i,1); %Optimal gains due to industrial policy
    %Now computing Harberger gains from trade policy
    lf_hat_opt = 1./(1-flab_real_opt(:,i));  %hat change in foreign labor when going from original equilibrium to optimum
    lf_hat_ip = 1./(1-flab_real_ip(:,i)); % hat change in foreign labor when going from original equilibrium to industrial policy only equilibrium
    flab_real_opt_ip=lf_hat_opt-lf_hat_ip;  % change in effective labor units supplied to foreigners when going from the industrial policy only equilibrium to optimum, relative to laissez-faire labor allocation
    flab_real_opt_ip(isnan(flab_real_opt_ip))=0;
    dw_harberger_trade(i,1)=(1/2)*sum(fs_initial(:,i).*flab_real_opt_ip.*-1./(theta+1)); %Harberger gains from optimal trade policy
    %Now computing Harberger gains from industrial policy
    lhat_opt = 1./(1-lab_real_opt(:,i));    %hat change in domestic labor allocation when going from original equilibrium to optimum
    lhat_tp=1./(1-lab_real_tp(:,i));    %hat change in domestic labor allocation when going from original equilibrium to equilibrium with trade policy only
    lab_real_opt_tp(:,i)=lhat_opt-lhat_tp;   %percent change in labor allocation when going from tp only equilibrium to optimum, relative to initial labor allocation
    dw_harberger_industrial(i,1)=(1/2)*sum(VA_initial(:,i).*lab_real_opt_tp(:,i).*gamma);
end

z=[countryid dw_harberger_opt dw_harberger_trade dw_harberger_industrial];
csvwrite('gains_soe_hte_ces_harberger.csv',z);

%Now saving initial share and reallocation terms to file for comparision
%for case without intermediates
m1 = [countryid VA_initial'];
csvwrite('va_initial.csv',m1);
m2 = [countryid lab_real_opt_tp'];
csvwrite('lab_real_opt_tp.csv',m2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from Industrial Policy at Autarky
%First compute hat welfare change when moving to autarky with no industrial
%policy
%Then compute hat change from moving to autarky with industrial policy
%Compare the two

%Het TE, Het SE (service 0 SE)
theta=h_theta; 
gamma = h_gamma;
counter=0;

for i=1:N
    %Defining data vectors and parameters
    V = V_i(i); %total value added
    X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
    D= D_i(i)/V; %Aggregate trade balance (divided by value added)
    beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
    lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
    V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
    
    VA_initial(:,i)=V_k; 
    fs_initial(:,i)=X;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Moving to Autarky with no trade policy
    %Specifying taxes and subsidies
    t_k=.95;
    s_k=zeros(K,1);
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    VA_TP=model_output.VA_share;
    gftp_hte_0(i,1)=model_output.welfare;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Now autarky with industrial policy
    %Specifying taxes and subsidies
    t_k=.95;
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gfip_hte_0(i,1)=real(model_output.welfare);
    counter=counter+1
end
gfoip_autarky = gfip_hte_0-gftp_hte_0

